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Abstract. - The velocity distribution of a granular gas is analyzed in terms of the Sonine 
polynomials expansion. We derive an analytical expression for the third Sonine coefficient 03. In 
contrast to frequently used assumptions this coefficient is of the same order of magnitude as the 
second Sonine coefficient a-z- For small inelasticity the theoretical result is in good agreement 
with numerical simulations. The next-order Sonine coefficients a^, a$ and a% are determined 
numerically. While these coefficients are negligible for small dissipation, their magnitude grows 
rapidly with increasing inelasticity for < e < 0.6. We conclude that this behavior of the 
Sonine coefficients manifests the break down of the Sonine polynomial expansion caused by the 
increasing impact of the overpopulated high-energy tail of the distribution function. 



Introduction. - The velocity distribution function of granular gases deviates from the 
Maxwell distribution, as first described by Goldshtein and Shapiro [1]. This deviation depends 
on the coefficient of restitution e, which quantifies the loss of energy for a collision of two 
particles i and j: 

1 + £ 1 + £ 

Vi=Vi — [(vi - vj) ■ e] e , = Vj + — — [(^ - Vj) ■ e] e. (1) 

Here v[ and stand for the post-collisional velocities where the unit vector of the relative 
particle position at the collision instant is e = (fi — fj) / \ fi — r)|. 

The deviation from the Maxwell distribution may be described by a Sonine polynomials 
expansion [1-3]. This expansion is applicable to the main part of the velocity distribution, 
excluding the high-energy tails, which is known to be exponentially overpopulated [4]. 

So far it was silently accepted that the Sonine expansion is a converging series while the 
exponential tail does not noticeably contribute to the coefficients of this expansion. This 
assumption was supported by Direct Simulation Monte Carlo (DSMC) of the Boltzmann 
equation [5], as well as by Molecular Dynamics simulations of Granular Gases [6]. Up to now, 

© EDP Sciences 



2 



EUROPHYSICS LETTERS 



however, neither the region of validity of the Sonine expansion is known, nor the impact of 
the exponential tail on the convergence of this series. 

Huthmann, Orza, and Brito [6] derived a system of equations for the Sonine coefficients 
from the Boltzmann equation and solved it perturbatively with the assumption that the Sonine 
coefficient at is of the order O (A fc ) with A being a small parameter. For e > 0.6 they obtained 
rapid decrease of the Sonine coefficients with escalating order k, while for smaller values of e 
the high-order coefficients, k > 3, were not negligible and could be of the same order as the 
first non-trivial coefficient a 2 - It was not evident, however, whether the Sonine polynomials 
expansion breaks down for e < 0.6, or the perturbative approach based on the conjecture 
A fe was inadequate. 

In the present study we derive the third Sonine coefficient, a 3 , and address the convergence 
of the Sonine polynomials expansion analytically and numerically by means of DSMC. We 
show that the Sonine coefficients do not decrease with increasing k > 3 for e < 0.6, i.e., the 
Sonine series diverges. We conclude that the breakdown of the Sonine expansion is caused by 
the increasing impact of the exponentially overpopulated tail for large dissipation. 

Sonine polynomials expansion. We consider a Granular Gas of particles of mass m 
which interact with a constant coefficient of restitution £=const. and neglect their rotational 
degrees of freedom. We assume that the gas is in the homogeneous cooling state at number 
density n — N/V. Moreover, we assume that the velocity distribution function of the gas 
f(v,t) has acquired its scaling form [1,4], i.e. 

where d is the system dimension and Vt{P) is the thermal velocity due to the temperature T(t), 

d „ f mv 2 , d mv 2 
2 

For e = const, the Boltzmann equation reduces to two uncoupled equations, for the tempera- 
ture and for the reduced distribution function /(c) [1,4]. The first equation reads 

dT 2 

-^ = -^92{<j)<J d ~ 1 nv T T i x 2l (4) 

where a is the particle diameter. For d = 3, the contact value of the pair distribution function 
is given by 32(0") = (1 — r//2) /(l — r/) 3 with the packing fraction i] = nn<T 3 /6 7 and 

lh = ~ j dcicfl{fj) (5) 
denotes the moments of the dimensionless collision integral [2,3,7], 



I m f mv 2 „, . d mv% 



I (/, /) = j dc 2 j deQ (-c 12 • e ) |c 12 • e I 



^/(cT)/(c2')-M)M) 



(6) 



The unit step function 0(a;) guarantees that only approaching particles collide, |ci2 • e | gives 
the length of the collision cylinder and c" and C2' denote the reduced velocities for the inverse 
collision, i.e., for the collision which results at the reduced velocities c\ and c*2. 
The second equation for the reduced distribution function reads 

1,2 ^ + c 1 ^-)f(c 1 )=I(f,f). (7) 



d V dci 
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For the case of elastic collisions, s = 1, the resulting velocity distribution is the Maxwell 
distribution. Therefore, for sufficiently large coefficient of restitution, e — > 1, we may assume 
that /(c) is close to the Maxwellian, 0(c) = 7r d / 2 exp (— c 2 ) . This suggests to expand the 
(unknown) distribution function /(c) around </>(c) in terms of orthogonal polynomials S p (x): 



/(c) = 0(cMc) 



1 + a,pS p (c 2 ) 
P=i 



(8) 



We chose S p (x) to be the Sonine polynomials (see, e.g., [8]): 



^ (n + l/2)!(p~n)!n! V j 

The first few of them, relevant for this study, read 
Si(x) = —x + —d 

S 2 (x) = ^x 2 -^(d + 2)x+^d(d + 2) (10) 
S 3 (x) = -^x 3 + J(d + A)x 2 - g(d + 2)(d + 4)1 + -^d(d + 2)(d + 4) . 

The expansion coefficients characterize the deviation of the distribution function from the 
Maxwell distribution, namely, they quantify the deviation of the moments of the distribution 
function, (c p ) = J dcc p f(c), from the corresponding values for the Maxwell distribution. The 
equations for the moments may be found multiplying both sides of Eq. (JJJ) by c\, integrating 
over ci, and using the orthogonality of the Sonine polynomials. This yields an infinite set of 
equations [2,3,8], 

dfi p = n 2 p (c p ) , p = 2,A,... (11) 

Since (c p ) and n P are expressed in terms of the Sonine coefficients, this set of equations can 
be used to determine the Sonine coefficients and, thus, to find the velocity distribution. To 
close the set of equations, a cutoff of the series (JSJ is applied, that is, it is assumed that the 
Sonine coefficients ak with k > fco are negligible. 

The first Sonine coefficient vanishes, a\ = 0, according to the definition of temperature [1]. 
Since (c 2 ) = the first equation in Eq. Ijllfl f° r P = 2 gives the identity. The second 
equation in Eq. (|llf) for p = 4 allows to find a<i by expressing (c 4 ), /i2 and /i4 in terms of 
ci2 and neglecting all other Sonine coefficients 03, 04, . . . [1]. Therefore, the first non-trivial 
Sonine coefficient is a 2 . 

Second and third Sonine coefficients. - We write /(c) = <p(c) [l + a2>5'2(c 2 ) + asS^c 2 )] 
for the velocity distribution function by assuming that 04, 05, . . . are negligible. Further, we 
write Eq. i|ll|) for p = 4 and p = 6 and express all quantities in these two equations in terms 
of <Z2 and 03. In particular, 

(c 4 ) = ~d(d + 2)(l + a 2 ), (c 6 ) = ld(d + 2){d + 4)(l + 3a 2 -a 3 ). (12) 

4 o 

To find the moments /12, /^4, and /ig we recast Eq. (0 using the properties of the collision 
integral (see, e.g., [8]): 

fj, p = -- I dc x I dc 2 I deQ {-c 12 ■ e) \c 12 ■ e | /(ci)/(c 2 ) A (c? + c%) , (13) 



4 



EUROPHYSICS LETTERS 



where A?p(ci) = ip (5j ) denotes the variation of some quantity ip (c) in a direct collision. 

The evaluation of integrals of the form of Eq. l|13fl has been described in detail in [8] . Using 
this approach, analytical expressions for /i 2l /J4, and [i§ may be obtained, which are rather 
cumbersome. Since it is expected that a 2 ^> a|, 03 3> a 2 , a 2 ^> 0203, and 03 ^> 0203, we keep 
in these moments only linear terms with respect to 02 and 03: 



T d/2 



Me 



27iT(d/2) 



27rr(d/2) 



r d/2 



27rr(d/2) 



^1 



1 - 

T 2 a2 -+ 
- D 2 a 2 



l a2 + k ai 



T303] 

h #303] 



where 



£>2 



(l- £ 
3 



2 +£ 



32 

_ (l-£ 2 ) 
128 
3 , 2 



d- 

(1 - e 2 ) (lOd + 10e 2 + 39) + (1 + e) (d - 1) 



10e 2 



97) _(l±4^i) (21 _ 5e) 



64 



(d- 



2e z 



19 

T 



(l-e 2 ) [l289-4(d + e 2 

256 v ' 1 v 



^3 



^1-e 2 ) [2537 + 4 (d + e ; 

2(d 2 



311 + 70,: 
! ) (583 



172d 2 



70s 2 ) + 236d 2 l B(e) 

16 



4e 2 ) 



(14) 



(15) 



1024 

(1 + e) [(d-3) (3- 

For 03 = 0, the above equations for \x 2 and /14 coincide with the equivalent equations in [2,3]. 

Substituting Eqs. (HJ, JUJ and into Eqs. (HJ for p = 4 and p = 6we obtain the 
Sonine cocfhcients a 2 and 03 in linear approximation. For d = 3, the result reads 

16 

a 2 = -— (240e 8 - 480e 7 + 3312e 6 - 7424e 5 + 3510e 4 - 364e 3 + 895e 2 + 1934e - 1623) 

(16) 



w 

128 



«3 



b(s) 



80e 8 - 160e 7 + 816e 6 - 1600e 5 + 154e 4 + 1548e 3 - 669e 2 - 386e + 217) 



6(e) = 2800£ 8 -5600e 7 +34800e 6 -84480£ 5 -4410e 4 +25716£ 3 +112155e 2 -172458e+214357. 

DSMC simulations. - To check the predictions of the theory and to study the behavior 
of higher Sonine coefficients we perform Monte Carlo simulations of the Boltzmann equation 
(DSMC) using 2 x 10 7 particles of unit mass. The coefficient of restitution was varied in the 
interval e G (0.1, 1) in steps of 0.01. To obtain smooth data, for each value of e we performed 
80 simulations and recorded the velocities of the particles when the system had reached a 
state with a scaling distribution function, Eq. J2J). From these snapshots we computed the 
temperature and the moments of the reduced distribution functions evaluating the averages 

1 N 1 JV / - N k 

r-5f2>«. <^>-sS(tBt) ■ (17) 

1 = 1 '4 = 1 V ' 
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Fig. 1 - The second Sonine coefficient over the coefficient of restitution e as given by Eq. JT^ 
where a-z and az are taken into account, 02(e) from the linear theory, where only 02 is taken into 
account [2], and 02(e) from the corresponding non- linear theory [3], together with DSMC results. 
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Fig. 2 - The coefficient 03 over the coefficient of restitution e due to Eq. Ijlrty and to DSMC. 
The Sonine coefficients can be computed using the relation (see e.g. [8]) 

Figure [2 shows a 2 as given by Eq. (|16() (which takes 0,3 into account), 0,2 as it follows 
from the linear theory [2] and CZ2 due to a non-linear theory [3] together with DSMC results. 
All approaches agree fairly well with the simulation results for small inelasticity, e < 0.6, and 
deviate noticeably for larger dissipation. Figure [5] shows 0,3 due to Eq. (|16|) together with the 
DSMC results. Again we see that the predictions of the new theory are in good agreement 
with the numerical results for e > 0.6. Similar numerical investigations of (13 have been done 
by Brey et al. [5] for e > 0.7 and of 03, CI4, a 5 by Nakanishi [9] for e > 0.9. 

Higher Sonine coefficients. In Fig. [3] we present the high-order Sonine coefficients 
as functions of the coefficient of restitution 04(e), 05(e), and ae(e). These coefficients are 
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Fig. 3 - High-order Sonine coefficients as functions of the coefficient of restitution (symbols). The 
lines guide the eye. 



very small for e < 0.6, which indicates the convergence of the Sonine expansion. For larger 
inelasticity, however, the absolute values of a p , p = 4, 5, 6 increase with increasing p. Hence 
we conclude that for large inelasticity, e < 0.6, the Sonine expansion does not converge. A 
similar result was reported in [6] , where the Sonine coefficients were found from the Boltzmann 
equation under the assumption au ^ \ k (Ac 1). 

What is the reason for the breakdown of the Sonine expansion with increasing inelasticity? 
This happens due to the increasing impact of the overpopulated tail of the velocity distribution 
function [4], which reads for c ^> 1 [2,4]: 

m - B ^ bc ; b = r(d+ i )n ■ ( 19 ) 

where \ii has been defined above, while the prefactor B is unknown. For small e the expo- 
nential overpopulation starts for velocities that are not significantly larger than the thermal 
velocity. In this case the contribution to the moments (c 2fe ) from the exponential tail rapidly 
grows with increasing k, which entails the corresponding growth of and ultimately, the 
breakdown of the Sonine expansion. Simulations show that at least in some range of e the 
main part of the distribution with c ~ 1 has a crossover to the tail distribution, Eq. I|19|l . at 
c ~ c* [10]. Assuming that the overpopulation of the tail starts at c* ~ b ~ 1/(1 — e 2 ) [8] one 
can analyze the contribution from the exponential tail to the moments (c 2fe ). 

The breakdown of the Sonine expansion may be also understood from a simple mathe- 
matical argument: The Sonine expansion, Eq. JSJ), contains only even powers of the scaled 
velocity c. The tail of the distribution decaying as exp(— c), c € (c*,oo), however, cannot be 
represented by a series in even powers of c. Therefore, for any value of e < 1, the presence of 
the asymptotic exponential tail must invalidate the Sonine expansion in the limit c — > oo. 

For e < 1 the tail starts at rather large velocity c* 3> 1, thus, only high-order terms of the 
Sonine expansion which are sensitive to large values of c are affected and the corresponding 
high-order Sonine coefficients a p are large. Indeed, for c> 1, according to Eq. JHJ 

<p(c) = ^-=J2 * P S P {c 2 ) - B exp (-be + c 2 )^B exp (c 2 ) = B £ ^ , (20) 
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while the Sonine polynomials may be approximated by their leading-order terms: 

(-1)p c 2 p 

S p «-^— ^ , thus a p ^(-l) p B for (21) 

The distribution function itself, however, is extremely small for these velocities. For larger 
inelasticity, that is, smaller e, the value of c* is not significantly larger than the thermal 
velocity. In this case, as shown in the present paper, already smaller Sonine coefficients such 
as <23 to ae may diverge. 

Nevertheless, the Sonine expansion remains a valuable tool for describing the main part of 
the velocity distribution, /(c), for c < c* which is the range of interest in most cases. When 
considering the expansion up to a®, the range e £ (0.6,1) seems to be a safe interval for 
applying the Sonine expansion. For higher-order expansions, however, the range of validity 
may be restricted to a smaller interval. A more quantitative discussion will be given in [10] 

Conclusion. We derived analytical expressions for the first two non-trivial Sonine 
coefficients as functions of the coefficient of restitution, 02(e) and 03(e), and show that the 
coefficient 03 is not negligible as compared to 02 as it was assumed in previous theories. We 
show that for small inelasticity, 0.6 < e < 1, au with k > 4 are significantly smaller than 
02, 03 and may be neglected, that is, the Sonine expansion converges. For this interval of e 
the obtained theoretical values of a 2 and a 3 are in a good agreement with numerical results. 
We also find numerically the coefficients 04, 05, and a 6 for a wide range of the restitution 
coefficient, 0.1 < e < 1. For large inelasticity, < e < 0.6, the high-order Sonine coefficient 
are of the same order as 02, a$ and, hence, may not be neglected. The reason for this behavior 
is the significant contribution of the exponentially overpopulated tail to the moments of the 
velocity distribution function. This contribution rapidly growths with increasing inelasticity 
and the order of the moments, ultimately undermining the Sonine polynomial expansion. 
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